Code
library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)Study 1
library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)Murphy (https://osf.io/3m5nh)
df1a <- haven::read_sav("../data/Murphy2020/Study 1.sav") |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))
df1b <- haven::read_sav("../data/Murphy2020/Study 6 IAS.sav") |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))Gaggero (2021) (https://osf.io/5x9sg)
load("../data/Gaggero2021/data.rda")
df2 <- data |>
mutate(Gender = as.character(gender)) |>
rename(
Age=age, Heart = `IAS 1`, Hungry = `IAS 2`, Breathing = `IAS 3`, Thirsty = `IAS 4`,
Urinate = `IAS 5`, Defecate = `IAS 6`, Taste = `IAS 7`, Vomit = `IAS 8`, Sneeze = `IAS 9`,
Cough = `IAS 10`, Temp = `IAS 11`, Sex_arousal = `IAS 12`, Wind = `IAS 13`, Burp = `IAS 14`,
Muscles = `IAS 15`, Bruise = `IAS 16`, Pain = `IAS 17`, Blood_Sugar = `IAS 18`,
Affective_touch = `IAS 19`, Tickle = `IAS 20`, Itch = `IAS 21`)
rm(data)Campos (2022) - Study 1 (https://osf.io/j6ef3)
df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |>
rename(
Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
)Todd (2022) - https://osf.io/ms354/
df4 <- haven::read_sav("../data/Todd2022/CompleteData_new.sav") |>
rename_with(\(x) str_remove(x, "IAS_"), .cols = starts_with("IAS_")) |>
rename(
Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
)Arslanova (2022) - https://osf.io/mp3cy/
df5 <- readxl::read_excel("../data/Arslanova2022/Murphy_Iacc_scale.xlsx") |>
filter(str_detect(Question.Key, "IAC_")) |>
filter(str_detect(Question.Key, "-quantised")) |>
pivot_wider(names_from = Question.Key, values_from = Response) |>
rename_with(\(x) str_remove(x, "IAC_"), .cols = starts_with("IAC_")) |>
rename_with(\(x) str_remove(x, "-quantised"), .cols = everything()) |>
rename(Sex_arousal = SexuallyAroused, Itch = Itchy, Sex_arousal = SexuallyAroused,
Temp = HotCold, Tickle = Ticklish, Breathing= BreatheFast,
Affective_touch= PleasantAffectionate, Blood_Sugar= BloodSugar,
Muscles=TiredMuscles, Heart= FastHeart, Taste=Tastes) |>
select(-starts_with("Participant")) |>
mutate_all(as.numeric)Brand (2022) - https://osf.io/xwz6g/
df6 <- haven::read_sav("../data/Brand2022/LatentVariableApproach.sav") |>
select(-SERIAL) |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |>
rename(
Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
)Brand (2023) - https://osf.io/3f2h6/
df7a <- haven::read_sav("../data/Brand2023/Data_Giessen.sav") |>
select(-COUNTRY_OTHER) |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(SEX == 1, "Male", ifelse(SEX == 2, "Female", "Other")))) |>
rename(
Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
)
df7b <- haven::read_sav("../data/Brand2023/Data_Mainz.sav") |>
select(-COUNTRY_OTHER, -EDUCATION_OTHER) |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |>
rename(
Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
)
df7c <- haven::read_sav("../data/Brand2023/Data_PotVie.sav") |>
select(-ID) |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |>
rename(
Age=AGE, Heart = IA02_01, Hungry = IA02_02, Breathing = IA02_03, Thirsty = IA02_04,
Urinate = IA02_05, Defecate = IA02_06, Taste = IA02_07, Vomit = IA02_08, Sneeze = IA02_09,
Cough = IA02_10, Temp = IA02_11, Sex_arousal = IA02_12, Wind = IA02_13, Burp = IA02_14,
Muscles = IA02_15, Bruise = IA02_16, Pain = IA02_17, Blood_Sugar = IA02_18,
Affective_touch = IA02_19, Tickle = IA02_20, Itch = IA02_21
)Lin (2023) - https://osf.io/3eztd - Note: No tickle because same Chinese character
df8a <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |>
select(-sex) |>
mutate_all(as.numeric) |>
mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |>
rename(
Age = age,
Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
Affective_touch = Touch, Itch = ITCH
)
df8b <- haven::read_sav("../data/Lin2023/Study 2.sav") |>
mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |>
rename(
Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
Affective_touch = Touch, Itch = ITCH)VonMohr (2023) - Study 3 (https://osf.io/7p9u5/)
df9 <- read.csv("../data/VonMohr2023/DataSet_study3.csv")
df9 <- df9[complete.cases(select(df9, starts_with("IAS_"))), ]
df9 <- df9 |>
mutate(Gender = as.character(ifelse(GENDER == "Man", "Male", ifelse(GENDER == "Woman", "Female", "Other")))) |>
rename(
Age=AGE, Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
)Makowski
df10a <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
rename(
Gender = Sex, Heart = Item_IAS_Interoception_1, Hungry = Item_IAS_Interoception_2,
Breathing = Item_IAS_Interoception_3, Thirsty = Item_IAS_Interoception_4,
Temp = Item_IAS_Interoception_5, Sex_arousal = Item_IAS_Interoception_6,
Urinate = Item_IAS_Elimination_1, Defecate = Item_IAS_Elimination_2,
Vomit = Item_IAS_Elimination_3, Wind = Item_IAS_Expulsion_1,
Burp = Item_IAS_Expulsion_2, Sneeze = Item_IAS_Expulsion_3,
Muscles = Item_IAS_Nociception_1, Bruise = Item_IAS_Nociception_2,
Pain = Item_IAS_Nociception_3, Affective_touch = Item_IAS_Skin_1,
Tickle = Item_IAS_Skin_2, Itch = Item_IAS_Skin_3
) |>
filter(!is.na(Urinate))
df10b <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
rename(
Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
)
data <- list(df1a=df1a, df1b=df1b, df2=df2, df3=df3, df4=df4, df5=df5, df6=df6, df7a=df7a, df7b=df7b, df7c=df7c, df8a=df8a, df8b=df8b, df9=df9, df10a=df10a, df10b=df10b)
save(data, file = "../data/data.RData")Total N = 31317.
The items with the differing distribution pattern (with a lower mode around 2/5) are “Affective Touch”, “Blood Sugar” and “Bruise”.
vars <- c(
"Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
"Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch"
)
dens1a <- estimate_density(select(df1a, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 1a")
dens1b <- estimate_density(select(df1b, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 1b")
dens2 <- estimate_density(select(df2, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 6")
dens7a <- estimate_density(select(df7a, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 7a")
dens7b <- estimate_density(select(df7b, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 7b")
dens7c <- estimate_density(select(df7c, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 7c")
dens8a <- estimate_density(select(df8a, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
mutate(Sample = "Sample 8a")
dens8b <- estimate_density(select(df8b, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
mutate(Sample = "Sample 8b")
dens9 <- estimate_density(select(df9, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 9")
dens10a <- estimate_density(select(df10a, all_of(setdiff(vars, c("Taste", "Cough", "Blood_Sugar")))), method = "kernSmooth") |>
mutate(Sample = "Sample 10a",
x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
dens10b <- estimate_density(select(df10b, all_of(vars)), method = "kernSmooth") |>
mutate(Sample = "Sample 10b",
x = datawizard::rescale(x, to=c(1, 5), range=c(0, 1)))
p1 <- rbind(dens1a, dens1b, dens2, dens3, dens4, dens5, dens6, dens7a, dens7b, dens7c, dens8a, dens8b, dens9, dens10a, dens10b) |>
mutate(Sample = fct_relevel(Sample, "Sample 1a", "Sample 1b", "Sample 2", "Sample 3", "Sample 4", "Sample 5", "Sample 6", "Sample 7a", "Sample 7b", "Sample 7c", "Sample 8a", "Sample 8b", "Sample 9", "Sample 10a", "Sample 10b")) |>
# mutate(Dashed = ifelse(Parameter %in% c("Affective_touch", "Blood_Sugar", "Bruise"), TRUE, FALSE)) |>
ggplot(aes(x = x, y = y, color = Parameter)) +
geom_line() +
scale_color_material_d() +
scale_linetype_manual(values = c("solid", "dashed"), guide="none") +
theme_minimal() +
theme(axis.text.y = element_blank(),
plot.title = element_text(face="bold")) +
labs(title = "Distribution of responses for all items across various datasets", x = "Response", y = "Density", color = "Item") +
facet_wrap(~Sample, scales = "free_y", nrow = 5)
ggsave("figures/Figure1.png", p1, width=10, height=10, dpi=200, bg="white")
p1data1a <- normalize(select(df1a, all_of(dens1a$Parameter)), range = c(1, 5))
data1b <- normalize(select(df1b, all_of(dens1b$Parameter)), range = c(1, 5))
data2 <- normalize(select(df2, all_of(dens2$Parameter)), range = c(1, 5))
data3 <- normalize(select(df3, all_of(dens3$Parameter)), range = c(1, 5))
data4 <- normalize(select(df4, all_of(dens4$Parameter)), range = c(1, 5))
data5 <- normalize(select(df5, all_of(dens5$Parameter)), range = c(1, 5))
data6 <- normalize(select(df6, all_of(dens6$Parameter)), range = c(1, 5))
data7a <- normalize(select(df7a, all_of(dens7a$Parameter)), range = c(1, 5))
data7b <- normalize(select(df7b, all_of(dens7b$Parameter)), range = c(1, 5))
data7c <- normalize(select(df7c, all_of(dens7c$Parameter)), range = c(1, 5))
data8a <- normalize(select(df8a, all_of(dens8a$Parameter)), range = c(1, 5))
data8b <- normalize(select(df8b, all_of(dens8b$Parameter)), range = c(1, 5))
data9 <- normalize(select(df9, all_of(dens9$Parameter)), range = c(1, 5))
data10a <- select(df10a, all_of(dens10a$Parameter))
data10b <- select(df10b, all_of(dens10b$Parameter))
data_all <- rbind(
data1a, data1b, data2, data3, data4, data5, data6, data7a, data7b, data7c,
mutate(data8a, Tickle = NA), mutate(data8b, Tickle = NA), data9,
mutate(data10a, Taste = NA, Cough = NA, Blood_Sugar = NA), data10b
) An overall positive intercorrelation pattern, with no clear structure emerging.
make_correlation <- function(df) {
correlation::correlation(df, redundant = TRUE) |>
correlation::cor_sort() |>
correlation::cor_lower() |>
mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
mutate(Parameter2 = fct_rev(Parameter2)) |>
ggplot(aes(x = Parameter1, y = Parameter2)) +
geom_tile(aes(fill = r), color = "white") +
geom_text(aes(label = val), size = 3) +
labs(title = "Correlation Matrix") +
scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks = FALSE)) +
theme_minimal() +
theme(
legend.title = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1)
)
}
make_correlation(data_all)make_correlation(data1a)make_correlation(data1b)make_correlation(data2)make_correlation(data3)make_correlation(data4)make_correlation(data5)make_correlation(data6)make_correlation(data7a)make_correlation(data7b)make_correlation(data7c)make_correlation(data8a)make_correlation(data8b)make_correlation(data9)make_correlation(data10a)make_correlation(data10b)Exploratory Graph Analysis (EGA) is a recently developed framework for psychometric assessment, that can be used to estimate the number of dimensions in questionnaire data using network estimation methods and community detection algorithms, and assess the stability of dimensions and items using bootstrapping. See https://r-ega.net/articles/ega.html for details.
Unique Variable Analysis (Christensen, Garrido, & Golino, 2023) uses the weighted topological overlap measure (Nowick et al., 2009) on an estimated network. Values greater than 0.25 are determined to have considerable local dependence (i.e., redundancy) that should be handled (variables with the highest maximum weighted topological overlap to all other variables (other than the one it is redundant with) should be removed).
uva <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Tickle Itch 0.368
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Wind Burp 0.293
Urinate Defecate 0.270
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Sneeze Cough 0.234
Heart Breathing 0.225
Hungry Thirsty 0.218
uva$keep_remove$keep
[1] "Itch"
$remove
[1] "Tickle"
uva <- EGAnet::UVA(data = data1a, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Tickle Itch 0.286
Wind Burp 0.275
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Sneeze Cough 0.244
Urinate Defecate 0.241
uva$keep_removeNULL
uva <- EGAnet::UVA(data = data1b, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Tickle Itch 0.350
Sneeze Cough 0.317
Wind Burp 0.309
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Urinate Defecate 0.278
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
uva$keep_remove$keep
[1] "Cough" "Burp" "Itch"
$remove
[1] "Sneeze" "Wind" "Tickle"
uva <- EGAnet::UVA(data = data3, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Tickle Itch 0.445
Sneeze Cough 0.318
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Wind Burp 0.256
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Bruise Blood_Sugar 0.219
Urinate Defecate 0.217
Heart Breathing 0.208
uva$keep_remove$keep
[1] "Sneeze" "Itch"
$remove
[1] "Cough" "Tickle"
uva <- EGAnet::UVA(data = data4, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Heart Breathing 0.415
Tickle Itch 0.351
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Hungry Thirsty 0.293
Urinate Defecate 0.282
Wind Burp 0.275
Sneeze Cough 0.251
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Muscles Pain 0.205
uva$keep_remove$keep
[1] "Heart" "Itch"
$remove
[1] "Breathing" "Tickle"
uva <- EGAnet::UVA(data = data5, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Tickle Itch 0.251
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Bruise Itch 0.241
uva$keep_removeNULL
uva <- EGAnet::UVA(data = data6, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Urinate Defecate 0.416
Sneeze Cough 0.332
Wind Burp 0.314
Tickle Itch 0.303
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Hungry Thirsty 0.274
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Heart Breathing 0.237
uva$keep_remove$keep
[1] "Urinate" "Sneeze" "Wind" "Itch"
$remove
[1] "Defecate" "Cough" "Burp" "Tickle"
uva <- EGAnet::UVA(data = data7a, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Sneeze Cough 0.434
Tickle Itch 0.376
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Wind Burp 0.295
Urinate Defecate 0.294
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
uva$keep_remove$keep
[1] "Sneeze" "Tickle"
$remove
[1] "Cough" "Itch"
uva <- EGAnet::UVA(data = data7b, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Urinate Defecate 0.317
Tickle Itch 0.311
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Heart Breathing 0.281
Sneeze Cough 0.261
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Wind Burp 0.237
uva$keep_remove$keep
[1] "Defecate" "Itch"
$remove
[1] "Urinate" "Tickle"
uva <- EGAnet::UVA(data = data7c, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Tickle Itch 0.457
Sneeze Cough 0.302
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Wind Burp 0.280
Urinate Defecate 0.254
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Hungry Thirsty 0.242
Vomit Sneeze 0.224
Heart Breathing 0.218
uva$keep_remove$keep
[1] "Cough" "Itch"
$remove
[1] "Sneeze" "Tickle"
uva <- EGAnet::UVA(data = data8a, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Urinate Defecate 0.237
Heart Breathing 0.218
Hungry Thirsty 0.213
uva$keep_removeNULL
uva <- EGAnet::UVA(data = data8b, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Urinate Defecate 0.277
Heart Breathing 0.267
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Wind Burp 0.206
uva$keep_removeNULL
uva <- EGAnet::UVA(data = data9, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
node_i node_j wto
Tickle Itch 0.379
Wind Burp 0.373
Urinate Defecate 0.341
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Heart Breathing 0.285
Sneeze Cough 0.278
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Hungry Thirsty 0.24
uva$keep_remove$keep
[1] "Defecate" "Wind" "Itch"
$remove
[1] "Urinate" "Burp" "Tickle"
uva <- EGAnet::UVA(data = data10a, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
node_i node_j wto
Tickle Itch 0.297
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
node_i node_j wto
Heart Breathing 0.209
uva$keep_removeNULL
uva <- EGAnet::UVA(data = data10b, cut.off = 0.3)
uvaVariable pairs with wTO > 0.30 (large-to-very large redundancy)
----
Variable pairs with wTO > 0.25 (moderate-to-large redundancy)
----
Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
uva$keep_removeNULL
make_egas <- function(data) {
egas <- list()
for (model in c("glasso")) { # , "TMFG"
for (algo in c("walktrap", "louvain")) {
for (type in c("ega", "ega.fit")) { # "" # "hierega"
if (type == "ega.fit" & algo == "louvain") next # Too slow
egas[[paste0(model, "_", algo, "_", type)]] <- EGAnet::bootEGA(
data = data,
seed = 123,
model = model,
algorithm = algo,
EGA.type = type,
type = "resampling",
plot.itemStability = FALSE,
verbose = FALSE
)
}
}
}
p <- EGAnet::compare.EGA.plots(
egas$glasso_walktrap_ega,
# egas$glasso_walktrap_ega.fit,
egas$glasso_louvain_ega,
# egas$TMFG_louvain_ega,
# egas$TMFG_walktrap_ega,
# egas$TMFG_walktrap_ega.fit,
labels = c(
"glasso_walktrap_ega",
# "glasso_walktrap_ega.fit",
"glasso_louvain_ega"
# "TMFG_louvain_ega",
# "TMFG_walktrap_ega",
# "TMFG_walktrap_ega.fit"
),
rows = 3,
plot.all = FALSE
)$all
list(egas = egas, p = p)
}
egas_all <- make_egas(data_all)
egas_all$pegas1a <- make_egas(data1a)
egas1a$pegas1b <- make_egas(data1b)
egas1b$pegas2 <- make_egas(data2)
egas2$pegas3 <- make_egas(data3)
egas3$pegas4 <- make_egas(data4)
egas4$pegas5 <- make_egas(data5)
egas5$pegas6 <- make_egas(data6)
egas6$pegas7a <- make_egas(data7a)
egas7a$pegas7b <- make_egas(data7b)
egas7b$pegas7c <- make_egas(data7c)
egas7c$pegas8a <- make_egas(data8a)
egas8a$pegas8b <- make_egas(data8b)
egas8b$pegas9 <- make_egas(data9)
egas9$pegas10a <- make_egas(data10a)
egas10a$pegas10b <- make_egas(data10b)
egas10b$ppatchwork::wrap_plots(lapply(egas_all$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas1a$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas1b$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas2$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas3$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas4$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas5$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas6$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas7a$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas7b$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas7c$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas8a$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas8b$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas9$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas10a$egas, plot), nrow = 3)patchwork::wrap_plots(lapply(egas10b$egas, plot), nrow = 3)ega <- egas_all$egas$glasso_louvain_ega$EGA
plot(ega)get_scores <- function(data, ega) {
EGAnet::net.scores(data, ega)$scores$std.scores |>
as.data.frame() |>
setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))
}
ega_scores_1a <- get_scores(data1a, ega)
ega_scores_1b <- get_scores(data1b, ega)
ega_scores_2 <- get_scores(data2, ega)
ega_scores_3 <- get_scores(data3, ega)
ega_scores_4 <- get_scores(data4, ega)
ega_scores_5 <- get_scores(data5, ega)
ega_scores_6 <- get_scores(data6, ega)
ega_scores_7a <- get_scores(data7a, ega)
ega_scores_7b <- get_scores(data7b, ega)
ega_scores_7c <- get_scores(data7c, ega)
ega_scores_8a <- get_scores(data8a, ega)
ega_scores_8b <- get_scores(data8b, ega)
ega_scores_9 <- get_scores(data9, ega)
ega_scores_10a <- get_scores(data10a, ega)
ega_scores_10b <- get_scores(data10b, ega)n <- parameters::n_factors(data_all, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | CNG | CNG |
| 3 | Scree (SE) | Scree_SE |
| 4 | beta | Multiple_regression |
| 4 | Optimal coordinates | Scree |
| 4 | Parallel analysis | Scree |
| 4 | Kaiser criterion | Scree |
n <- parameters::n_factors(data1a, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Bentler | Bentler |
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | CNG | CNG |
| 3 | Scree (SE) | Scree_SE |
| 3 | BIC | BIC |
| 4 | beta | Multiple_regression |
| 4 | Optimal coordinates | Scree |
| 4 | Parallel analysis | Scree |
| 4 | Kaiser criterion | Scree |
n <- parameters::n_factors(data1b, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Bentler | Bentler |
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 3 | CNG | CNG |
| 3 | Scree (SE) | Scree_SE |
| 3 | BIC | BIC |
| 4 | beta | Multiple_regression |
| 4 | Optimal coordinates | Scree |
| 4 | Parallel analysis | Scree |
| 4 | Kaiser criterion | Scree |
| 5 | BIC (adjusted) | BIC |
n <- parameters::n_factors(data2, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | CNG | CNG |
| 4 | beta | Multiple_regression |
| 4 | Scree (SE) | Scree_SE |
| 4 | BIC | BIC |
| 5 | Optimal coordinates | Scree |
| 5 | Parallel analysis | Scree |
| 5 | Kaiser criterion | Scree |
n <- parameters::n_factors(data3, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 3 | CNG | CNG |
| 3 | Velicer’s MAP | Velicers_MAP |
| 4 | beta | Multiple_regression |
| 4 | Optimal coordinates | Scree |
| 4 | Parallel analysis | Scree |
| 4 | Kaiser criterion | Scree |
| 4 | Scree (SE) | Scree_SE |
| 4 | BIC | BIC |
n <- parameters::n_factors(data4, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | Bentler | Bentler |
| 3 | CNG | CNG |
| 4 | beta | Multiple_regression |
| 4 | Optimal coordinates | Scree |
| 4 | Parallel analysis | Scree |
| 4 | Kaiser criterion | Scree |
| 4 | Scree (SE) | Scree_SE |
| 5 | BIC | BIC |
n <- parameters::n_factors(data5, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Optimal coordinates | Scree |
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 1 | BIC | BIC |
| 3 | CNG | CNG |
| 3 | BIC (adjusted) | BIC |
| 4 | Bartlett | Barlett |
| 4 | beta | Multiple_regression |
| 4 | Scree (SE) | Scree_SE |
| 5 | Anderson | Barlett |
| 5 | Lawley | Barlett |
n <- parameters::n_factors(data6, n_max = 12)Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
ultra-Heywood case was detected. Examine the results carefully
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | Bentler | Bentler |
| 3 | CNG | CNG |
| 3 | Scree (SE) | Scree_SE |
| 4 | beta | Multiple_regression |
| 4 | Optimal coordinates | Scree |
| 4 | Parallel analysis | Scree |
| 4 | Kaiser criterion | Scree |
n <- parameters::n_factors(data7a, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Bentler | Bentler |
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | CNG | CNG |
| 3 | Scree (SE) | Scree_SE |
| 3 | BIC | BIC |
| 4 | beta | Multiple_regression |
| 4 | Optimal coordinates | Scree |
| 4 | Parallel analysis | Scree |
| 4 | Kaiser criterion | Scree |
n <- parameters::n_factors(data7b, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | CNG | CNG |
| 4 | beta | Multiple_regression |
| 4 | Scree (SE) | Scree_SE |
| 5 | Bentler | Bentler |
| 5 | Optimal coordinates | Scree |
| 5 | Parallel analysis | Scree |
| 5 | Kaiser criterion | Scree |
n <- parameters::n_factors(data7c, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | Bentler | Bentler |
| 3 | CNG | CNG |
| 4 | beta | Multiple_regression |
| 4 | Optimal coordinates | Scree |
| 4 | Scree (SE) | Scree_SE |
n <- parameters::n_factors(data8a, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | Bentler | Bentler |
| 3 | CNG | CNG |
| 4 | beta | Multiple_regression |
| 4 | Scree (SE) | Scree_SE |
| 5 | Optimal coordinates | Scree |
| 5 | Parallel analysis | Scree |
| 5 | Kaiser criterion | Scree |
| 5 | BIC | BIC |
n <- parameters::n_factors(data8b, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | CNG | CNG |
| 4 | beta | Multiple_regression |
| 4 | BIC | BIC |
| 5 | Optimal coordinates | Scree |
| 5 | Parallel analysis | Scree |
| 5 | Kaiser criterion | Scree |
| 5 | Scree (SE) | Scree_SE |
| 5 | BIC (adjusted) | BIC |
n <- parameters::n_factors(data9, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | CNG | CNG |
| 3 | Scree (SE) | Scree_SE |
| 4 | beta | Multiple_regression |
| 4 | Optimal coordinates | Scree |
| 4 | Parallel analysis | Scree |
| 4 | Kaiser criterion | Scree |
n <- parameters::n_factors(data10a, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | CNG | CNG |
| 3 | Scree (SE) | Scree_SE |
| 4 | beta | Multiple_regression |
| 4 | Optimal coordinates | Scree |
| 4 | Parallel analysis | Scree |
| 4 | Kaiser criterion | Scree |
| 4 | BIC | BIC |
| 4 | BIC (adjusted) | BIC |
n <- parameters::n_factors(data10b, n_max = 12)
plot(n)knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 3, 4, 5), ], row.names = FALSE)| n_Factors | Method | Family |
|---|---|---|
| 1 | Bentler | Bentler |
| 1 | Acceleration factor | Scree |
| 1 | Scree (R2) | Scree_SE |
| 1 | VSS complexity 1 | VSS |
| 1 | Velicer’s MAP | Velicers_MAP |
| 3 | CNG | CNG |
| 4 | beta | Multiple_regression |
| 4 | Optimal coordinates | Scree |
| 4 | Parallel analysis | Scree |
| 4 | Kaiser criterion | Scree |
| 4 | Scree (SE) | Scree_SE |
| 5 | BIC (adjusted) | BIC |
efa5 <- parameters::factor_analysis(data_all, n = 4, rotation = "oblimin", sort = TRUE)Loading required namespace: GPArotation
plot(efa5)display(efa5)| Variable | MR1 | MR3 | MR2 | MR4 | Complexity | Uniqueness |
|---|---|---|---|---|---|---|
| Burp | 0.74 | 0.03 | -0.02 | -0.07 | 1.02 | 0.48 |
| Cough | 0.70 | -0.02 | 0.02 | 0.06 | 1.02 | 0.47 |
| Wind | 0.67 | 0.02 | -0.02 | -2.95e-03 | 1.00 | 0.55 |
| Sneeze | 0.63 | -0.06 | 0.04 | 0.14 | 1.12 | 0.52 |
| Vomit | 0.44 | 0.06 | 0.04 | 0.18 | 1.39 | 0.63 |
| Temp | 0.27 | 0.24 | 0.05 | 0.22 | 2.98 | 0.61 |
| Sex_arousal | 0.26 | 0.22 | 0.03 | 0.20 | 2.90 | 0.67 |
| Taste | 0.24 | 0.20 | 0.10 | 0.15 | 2.99 | 0.70 |
| Breathing | 0.05 | 0.61 | -0.04 | 0.04 | 1.03 | 0.59 |
| Hungry | -0.11 | 0.54 | -5.62e-03 | 0.19 | 1.34 | 0.64 |
| Heart | 0.08 | 0.49 | 0.03 | -0.07 | 1.10 | 0.73 |
| Thirsty | -0.07 | 0.47 | -7.04e-03 | 0.25 | 1.60 | 0.64 |
| Pain | 0.16 | 0.40 | 0.10 | 0.10 | 1.60 | 0.62 |
| Muscles | 0.30 | 0.36 | 0.07 | 0.03 | 2.03 | 0.59 |
| Blood_Sugar | 0.11 | 0.35 | 0.22 | -0.18 | 2.47 | 0.76 |
| Bruise | 0.25 | 0.34 | 0.22 | -0.19 | 3.30 | 0.65 |
| Tickle | -0.04 | -0.04 | 0.82 | 0.04 | 1.02 | 0.37 |
| Itch | 0.02 | 6.24e-03 | 0.75 | -0.04 | 1.01 | 0.43 |
| Affective_touch | 0.08 | 0.19 | 0.33 | 0.11 | 2.00 | 0.70 |
| Urinate | 0.04 | 0.10 | 0.04 | 0.65 | 1.06 | 0.45 |
| Defecate | 0.15 | 0.04 | 0.05 | 0.61 | 1.14 | 0.47 |
The 4 latent factors (oblimin rotation) accounted for 41.65% of the total variance of the original data (MR1 = 14.63%, MR3 = 11.59%, MR2 = 8.13%, MR4 = 7.30%).
efa4 <- parameters::factor_analysis(data1a, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)display(efa4)| Variable | MR2 | MR3 | MR1 | MR4 | Complexity | Uniqueness |
|---|---|---|---|---|---|---|
| Tickle | 0.72 | 0.03 | -0.10 | 0.11 | 1.09 | 0.47 |
| Itch | 0.67 | 0.14 | -9.44e-03 | -0.10 | 1.13 | 0.46 |
| Affective_touch | 0.49 | -0.13 | 0.29 | 0.13 | 1.97 | 0.59 |
| Blood_Sugar | 0.29 | 0.21 | 0.18 | -0.23 | 3.51 | 0.75 |
| Burp | -0.03 | 0.73 | -0.03 | -1.21e-03 | 1.00 | 0.50 |
| Cough | 7.56e-03 | 0.69 | 7.00e-03 | -0.01 | 1.00 | 0.52 |
| Sneeze | 0.07 | 0.54 | 7.09e-03 | 0.20 | 1.30 | 0.55 |
| Wind | 0.08 | 0.46 | 0.07 | 0.16 | 1.37 | 0.62 |
| Hungry | -0.16 | 0.04 | 0.65 | 0.03 | 1.13 | 0.60 |
| Thirsty | -4.84e-03 | 0.01 | 0.49 | 0.11 | 1.11 | 0.70 |
| Breathing | 0.14 | 0.05 | 0.44 | 0.03 | 1.24 | 0.70 |
| Heart | 0.11 | -0.07 | 0.40 | 0.08 | 1.32 | 0.78 |
| Bruise | 0.31 | 0.25 | 0.34 | -0.25 | 3.74 | 0.59 |
| Pain | 0.29 | 0.01 | 0.32 | 0.20 | 2.67 | 0.61 |
| Temp | 0.19 | 0.13 | 0.29 | 0.16 | 2.91 | 0.67 |
| Muscles | 0.20 | 0.13 | 0.26 | 0.16 | 3.24 | 0.69 |
| Sex_arousal | 0.18 | 0.02 | 0.21 | 0.15 | 2.81 | 0.82 |
| Defecate | 0.10 | 0.11 | 2.77e-03 | 0.66 | 1.10 | 0.46 |
| Urinate | -0.11 | 0.13 | 0.24 | 0.52 | 1.65 | 0.54 |
| Taste | 0.23 | 0.06 | 0.17 | 0.30 | 2.60 | 0.69 |
| Vomit | 0.21 | 0.18 | 0.11 | 0.25 | 3.26 | 0.71 |
The 4 latent factors (oblimin rotation) accounted for 37.94% of the total variance of the original data (MR2 = 10.46%, MR3 = 10.17%, MR1 = 10.08%, MR4 = 7.23%).
efa4 <- parameters::factor_analysis(data1b, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)display(efa4)| Variable | MR1 | MR2 | MR3 | MR4 | Complexity | Uniqueness |
|---|---|---|---|---|---|---|
| Urinate | 0.69 | -0.18 | 0.13 | -0.04 | 1.22 | 0.54 |
| Defecate | 0.68 | -0.13 | 0.04 | 0.06 | 1.10 | 0.55 |
| Breathing | 0.56 | 0.03 | 0.06 | 7.96e-03 | 1.03 | 0.63 |
| Hungry | 0.55 | 0.15 | -7.15e-03 | 0.02 | 1.16 | 0.60 |
| Thirsty | 0.52 | 0.05 | 0.14 | -0.01 | 1.16 | 0.62 |
| Sex_arousal | 0.50 | 0.14 | -0.04 | 0.20 | 1.49 | 0.56 |
| Temp | 0.46 | 0.23 | 0.16 | -0.01 | 1.77 | 0.52 |
| Pain | 0.46 | 0.43 | -0.05 | -0.03 | 2.03 | 0.49 |
| Taste | 0.40 | 0.21 | -0.07 | 0.12 | 1.77 | 0.69 |
| Heart | 0.40 | 0.06 | -0.04 | 0.15 | 1.36 | 0.76 |
| Muscles | 0.35 | 0.24 | 0.07 | 0.09 | 2.07 | 0.65 |
| Vomit | 0.32 | 0.06 | 0.23 | 0.11 | 2.20 | 0.68 |
| Tickle | 0.02 | 0.69 | 0.01 | 0.04 | 1.01 | 0.47 |
| Itch | -0.02 | 0.67 | 0.18 | 0.03 | 1.16 | 0.41 |
| Bruise | 0.07 | 0.44 | 0.15 | 0.06 | 1.32 | 0.65 |
| Blood_Sugar | -0.08 | 0.38 | 0.16 | 0.09 | 1.61 | 0.77 |
| Affective_touch | 0.38 | 0.38 | -0.13 | 0.05 | 2.26 | 0.63 |
| Sneeze | 0.05 | 0.05 | 0.71 | 0.01 | 1.02 | 0.41 |
| Cough | 0.04 | 0.02 | 0.71 | 0.04 | 1.01 | 0.43 |
| Wind | 7.52e-03 | -0.04 | -5.18e-03 | 0.94 | 1.00 | 0.15 |
| Burp | -0.01 | 0.27 | 0.17 | 0.47 | 1.91 | 0.48 |
The 4 latent factors (oblimin rotation) accounted for 44.25% of the total variance of the original data (MR1 = 17.71%, MR2 = 11.72%, MR3 = 7.65%, MR4 = 7.17%).
We discarded the following items: - Tickle: not present in the same dataset and consistently flagged as redundant.
Ambiguous: - Temp - Vomit - Affective_touch - Sex_arousal - Taste
m1 <- "
Interoception =~ Hungry + Thirsty + Urinate + Defecate + Itch + Bruise + Muscles + Pain + Breathing + Heart + Cough + Sneeze + Wind + Burp
"
m4 <- "
Sustenance =~ Hungry + Thirsty + Urinate + Defecate + Muscles + Pain
Skin =~ Itch + Bruise
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"
m5 <- "
Sustenance =~ Hungry + Thirsty + Pain
UrinateDefecate =~ Urinate + Defecate
Skin =~ Itch + Bruise + Muscles
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"
m6 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"
m7 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
CoughSneeze =~ Cough + Sneeze
WindBurp =~ Wind + Burp
"
is_converged <- function(m) {
ev <- eigen(lavaan::lavTech(m, "cov.lv")[[1]], symmetric=TRUE, only.values=TRUE)$values
if (any(ev < -1 * .Machine$double.eps^(3/4))) {
"Not converged"
} else {
""
}
}
compare_mods <- function(data) {
m1 <- lavaan::cfa(m1, data = data)
m4 <- lavaan::cfa(m4, data = data)
m5 <- lavaan::cfa(m5, data = data)
m6 <- lavaan::cfa(m6, data = data)
m7 <- lavaan::cfa(m7, data = data)
anova(m1, m4, m5, m6, m7) |>
parameters::parameters() |>
mutate(AIC = format(round(AIC, 0), nsmall = 0),
BIC = format(round(BIC, 0), nsmall = 0),
Chi2 = format(round(Chi2, 2), nsmall = 2),
Converged = sapply(list(m1, m4, m5, m6, m7), is_converged))
}
display(compare_mods(data_all))| Parameter | AIC | BIC | Chi2 | df | p | Converged |
|---|---|---|---|---|---|---|
| m7 | -168983 | -168574 | 2080.93 | 56 | ||
| m6 | -164335 | -163977 | 6740.02 | 62 | < .001 | |
| m5 | -160047 | -159731 | 11038.16 | 67 | < .001 | |
| m4 | -154973 | -154690 | 16120.58 | 71 | < .001 | |
| m1 | -141854 | -141621 | 29251.50 | 77 | < .001 |
m7_all <- lavaan::cfa(m7, data = data_all)
parameters::parameters(m7_all, standardize=TRUE) |>
display()| Link | Coefficient | SE | 95% CI | z | p |
|---|---|---|---|---|---|
| HungryThirsty =~ Hungry | 0.68 | 4.79e-03 | (0.67, 0.69) | 142.07 | < .001 |
| HungryThirsty =~ Thirsty | 0.71 | 4.77e-03 | (0.70, 0.72) | 148.26 | < .001 |
| UrinateDefecate =~ Urinate | 0.77 | 3.94e-03 | (0.76, 0.78) | 195.43 | < .001 |
| UrinateDefecate =~ Defecate | 0.77 | 3.93e-03 | (0.77, 0.78) | 196.53 | < .001 |
| ItchBruise =~ Itch | 0.50 | 5.73e-03 | (0.49, 0.51) | 87.60 | < .001 |
| ItchBruise =~ Bruise | 0.67 | 5.99e-03 | (0.66, 0.68) | 112.01 | < .001 |
| MusclesPain =~ Muscles | 0.71 | 4.15e-03 | (0.70, 0.72) | 171.30 | < .001 |
| MusclesPain =~ Pain | 0.68 | 4.22e-03 | (0.68, 0.69) | 161.84 | < .001 |
| HeartBreathing =~ Breathing | 0.79 | 5.19e-03 | (0.78, 0.80) | 151.51 | < .001 |
| HeartBreathing =~ Heart | 0.59 | 5.13e-03 | (0.58, 0.60) | 114.83 | < .001 |
| CoughSneeze =~ Cough | 0.81 | 3.49e-03 | (0.80, 0.81) | 230.73 | < .001 |
| CoughSneeze =~ Sneeze | 0.75 | 3.65e-03 | (0.74, 0.76) | 205.71 | < .001 |
| WindBurp =~ Wind | 0.77 | 3.58e-03 | (0.76, 0.78) | 215.22 | < .001 |
| WindBurp =~ Burp | 0.82 | 3.44e-03 | (0.81, 0.83) | 238.24 | < .001 |
| Link | Coefficient | SE | 95% CI | z | p |
|---|---|---|---|---|---|
| HungryThirsty ~~ UrinateDefecate | 0.65 | 6.26e-03 | (0.63, 0.66) | 103.09 | < .001 |
| HungryThirsty ~~ ItchBruise | 0.45 | 9.03e-03 | (0.43, 0.46) | 49.41 | < .001 |
| HungryThirsty ~~ MusclesPain | 0.63 | 7.04e-03 | (0.61, 0.64) | 89.04 | < .001 |
| HungryThirsty ~~ HeartBreathing | 0.64 | 7.01e-03 | (0.63, 0.66) | 91.55 | < .001 |
| HungryThirsty ~~ CoughSneeze | 0.45 | 7.14e-03 | (0.44, 0.47) | 63.70 | < .001 |
| HungryThirsty ~~ WindBurp | 0.43 | 7.12e-03 | (0.42, 0.44) | 60.47 | < .001 |
| UrinateDefecate ~~ ItchBruise | 0.43 | 8.38e-03 | (0.42, 0.45) | 51.81 | < .001 |
| UrinateDefecate ~~ MusclesPain | 0.64 | 6.23e-03 | (0.63, 0.66) | 103.29 | < .001 |
| UrinateDefecate ~~ HeartBreathing | 0.53 | 6.81e-03 | (0.52, 0.55) | 78.24 | < .001 |
| UrinateDefecate ~~ CoughSneeze | 0.56 | 6.03e-03 | (0.54, 0.57) | 92.16 | < .001 |
| UrinateDefecate ~~ WindBurp | 0.50 | 6.19e-03 | (0.49, 0.52) | 81.27 | < .001 |
| ItchBruise ~~ MusclesPain | 0.80 | 7.96e-03 | (0.79, 0.82) | 100.88 | < .001 |
| ItchBruise ~~ HeartBreathing | 0.51 | 8.74e-03 | (0.49, 0.53) | 58.21 | < .001 |
| ItchBruise ~~ CoughSneeze | 0.64 | 7.69e-03 | (0.62, 0.65) | 82.74 | < .001 |
| ItchBruise ~~ WindBurp | 0.64 | 7.55e-03 | (0.62, 0.65) | 84.63 | < .001 |
| MusclesPain ~~ HeartBreathing | 0.64 | 6.99e-03 | (0.62, 0.65) | 91.42 | < .001 |
| MusclesPain ~~ CoughSneeze | 0.66 | 6.06e-03 | (0.65, 0.68) | 109.42 | < .001 |
| MusclesPain ~~ WindBurp | 0.63 | 6.10e-03 | (0.62, 0.65) | 103.72 | < .001 |
| HeartBreathing ~~ CoughSneeze | 0.52 | 6.81e-03 | (0.50, 0.53) | 75.70 | < .001 |
| HeartBreathing ~~ WindBurp | 0.45 | 6.95e-03 | (0.44, 0.47) | 65.21 | < .001 |
| CoughSneeze ~~ WindBurp | 0.72 | 4.84e-03 | (0.71, 0.73) | 149.24 | < .001 |
display(compare_mods(data1a))| Parameter | AIC | BIC | Chi2 | df | p | Converged |
|---|---|---|---|---|---|---|
| m7 | -1951 | -1749 | 88.22 | 56 | ||
| m6 | -1925 | -1748 | 126.29 | 62 | < .001 | |
| m5 | -1860 | -1703 | 201.38 | 67 | < .001 | Not converged |
| m4 | -1834 | -1694 | 234.82 | 71 | < .001 | |
| m1 | -1661 | -1546 | 420.24 | 77 | < .001 |
Conclusion: The 7th factor solution was preferred in most datasets, including by indices that penalize increased number of parameters (such as BIC).
m7h1 <- paste0(m7, "
Interoception =~ HungryThirsty + UrinateDefecate + ItchBruise + MusclesPain + HeartBreathing + CoughSneeze + WindBurp
")
m7h2 <- paste0(m7, "
Expulsion =~ CoughSneeze + WindBurp
Interoception =~ HungryThirsty + UrinateDefecate + HeartBreathing
")
m7h3 <- paste0(m7, "
Valenced =~ MusclesPain + HeartBreathing + ItchBruise
Expulsion =~ CoughSneeze + WindBurp
Homeostasis =~ HungryThirsty + UrinateDefecate
")
compare_mods_h <- function(data) {
m7 <- lavaan::cfa(m7, data = data)
m7h2 <- lavaan::cfa(m7h2, data = data)
m7h1 <- lavaan::cfa(m7h1, data = data)
m7h3 <- lavaan::cfa(m7h3, data = data)
anova(m7, m7h1, m7h2, m7h3) |>
parameters::parameters() |>
mutate(AIC = format(round(AIC, 0), nsmall = 0),
BIC = format(round(BIC, 0), nsmall = 0),
Chi2 = format(round(Chi2, 2), nsmall = 2),
Converged = sapply(list(m7, m7h1, m7h2, m7h3), is_converged))
}
display(compare_mods_h(data_all))| Parameter | AIC | BIC | Chi2 | df | p | Converged |
|---|---|---|---|---|---|---|
| m7 | -168983 | -168574 | 2080.93 | 56 | ||
| m7h2 | -168224 | -167899 | 2859.83 | 66 | < .001 | |
| m7h3 | -166932 | -166615 | 4153.51 | 67 | < .001 | |
| m7h1 | -164501 | -164210 | 6590.01 | 70 | < .001 |
No evidence for higher order factors.
display(compare_mods_h(data1a))| Parameter | AIC | BIC | Chi2 | df | p | Converged |
|---|---|---|---|---|---|---|
| m7 | -1951 | -1749 | 88.22 | 56 | ||
| m7h3 | -1904 | -1747 | 157.51 | 67 | < .001 | |
| m7h1 | -1873 | -1729 | 193.71 | 70 | < .001 |
rbind(
performance::performance(m7_all, metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "All"),
performance::performance(update(m7_all, data = data1a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 1a"),
performance::performance(update(m7_all, data = data1b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 1b"),
performance::performance(update(m7_all, data = data2), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 2"),
performance::performance(update(m7_all, data = data3), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 3"),
performance::performance(update(m7_all, data = data4), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 4"),
performance::performance(update(m7_all, data = data5), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 5"),
performance::performance(update(m7_all, data = data6), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 6"),
performance::performance(update(m7_all, data = data7a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 7a"),
performance::performance(update(m7_all, data = data7b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 7b"),
performance::performance(update(m7_all, data = data7c), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 7c"),
performance::performance(update(m7_all, data = data8a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 8a"),
performance::performance(update(m7_all, data = data8b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 8b"),
performance::performance(update(m7_all, data = data9), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 9"),
performance::performance(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data10a), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 10a"),
performance::performance(update(m7_all, data = data10b), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
mutate(Sample = "Sample 10b")
) |>
display()| Chi2 | RMSEA | CFI | SRMR | Sample |
|---|---|---|---|---|
| 2080.93 | 0.03 | 0.98 | 0.02 | All |
| 88.22 | 0.04 | 0.98 | 0.03 | Sample 1a |
| 66.50 | 0.02 | 0.99 | 0.03 | Sample 1b |
| 123.77 | 0.04 | 0.97 | 0.03 | Sample 2 |
| 114.03 | 0.04 | 0.98 | 0.03 | Sample 3 |
| 81.35 | 0.02 | 0.99 | 0.02 | Sample 4 |
| 83.52 | 0.05 | 0.92 | 0.06 | Sample 5 |
| 152.88 | 0.05 | 0.97 | 0.03 | Sample 6 |
| 106.90 | 0.04 | 0.97 | 0.04 | Sample 7a |
| 188.21 | 0.03 | 0.98 | 0.02 | Sample 7b |
| 122.56 | 0.04 | 0.97 | 0.03 | Sample 7c |
| 115.89 | 0.03 | 0.98 | 0.02 | Sample 8a |
| 83.61 | 0.03 | 0.98 | 0.03 | Sample 8b |
| 2003.31 | 0.04 | 0.98 | 0.02 | Sample 9 |
| 82.59 | 0.05 | 0.97 | 0.03 | Sample 10a |
| 91.47 | 0.05 | 0.97 | 0.04 | Sample 10b |
Note the usage of descriptive factor names relating directly to the items to avoid abstraction.
# Add empirical variables
add_empirical <- function(data, df) {
x <- data |>
mutate(
HungryThirsty = (df$Hungry + df$Thirsty) / 2,
UrinateDefecate = (df$Urinate + df$Defecate) / 2,
MusclesPain = (df$Muscles + df$Pain) / 2,
HeartBreathing = (df$Heart + df$Breathing) / 2,
ItchBruise = (df$Itch + df$Bruise) / 2,
WindBurp = (df$Wind + df$Burp) / 2)
if("Cough" %in% names(df)) {
mutate(x, CoughSneeze = (df$Cough + df$Sneeze) / 2)
} else {
mutate(x, CoughSneeze = df$Sneeze)
}
}
# Deal with missing data
temp <- predict(m7_all, newdata=data7c)
scores7c <- data.frame(matrix(ncol = 7, nrow = nrow(ega_scores_7c)))
scores7c[!is.na(ega_scores_7c$EGA1), ] <- temp
names(scores7c) <- names(as.data.frame(temp))
scores <- list(
sample1a = cbind(ega_scores_1a, data_addprefix(as.data.frame(predict(m7_all, newdata=data1a)), "CFA_")) |>
add_empirical(data1a) |>
mutate(Sample = "Sample 1a"),
sample1b = cbind(ega_scores_1b, data_addprefix(as.data.frame(predict(m7_all, newdata=data1b)), "CFA_")) |>
add_empirical(data1b) |>
mutate(Sample = "Sample 1b"),
sample2 = cbind(ega_scores_2, data_addprefix(as.data.frame(predict(m7_all, newdata=data2)), "CFA_")) |>
add_empirical(data2) |>
mutate(Sample = "Sample 2"),
sample3 = cbind(ega_scores_3, data_addprefix(as.data.frame(predict(m7_all, newdata=data3)), "CFA_")) |>
add_empirical(data3) |>
mutate(Sample = "Sample 3"),
sample4 = cbind(ega_scores_4, data_addprefix(as.data.frame(predict(m7_all, newdata=data4)), "CFA_")) |>
add_empirical(data4) |>
mutate(Sample = "Sample 4"),
sample5 = cbind(ega_scores_5, data_addprefix(as.data.frame(predict(m7_all, newdata=data5)), "CFA_")) |>
add_empirical(data5) |>
mutate(Sample = "Sample 5"),
sample6 = cbind(ega_scores_6, data_addprefix(as.data.frame(predict(m7_all, newdata=data6)), "CFA_")) |>
add_empirical(data6) |>
mutate(Sample = "Sample 6"),
sample7a = cbind(ega_scores_7a, data_addprefix(as.data.frame(predict(m7_all, newdata=data7a)), "CFA_")) |>
add_empirical(data7a) |>
mutate(Sample = "Sample 7a"),
sample7b = cbind(ega_scores_7b, data_addprefix(as.data.frame(predict(m7_all, newdata=data7b)), "CFA_")) |>
add_empirical(data7b) |>
mutate(Sample = "Sample 7b"),
sample7c = cbind(ega_scores_7c, data_addprefix(scores7c, "CFA_")) |>
add_empirical(data7c) |>
mutate(Sample = "Sample 7c"),
sample8a = cbind(ega_scores_8a, data_addprefix(as.data.frame(predict(m7_all, newdata=data8a)), "CFA_")) |>
add_empirical(data8a) |>
mutate(Sample = "Sample 8a"),
sample8b = cbind(ega_scores_8b, data_addprefix(as.data.frame(predict(m7_all, newdata=data8b)), "CFA_")) |>
add_empirical(data8b) |>
mutate(Sample = "Sample 8b"),
sample9 = cbind(ega_scores_9, data_addprefix(as.data.frame(predict(m7_all, newdata=data9)), "CFA_")) |>
add_empirical(data9) |>
mutate(Sample = "Sample 9"),
sample10a = cbind(ega_scores_10a, data_addprefix(as.data.frame(predict(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data_all), newdata=data10a)), "CFA_")) |>
add_empirical(data10a) |>
mutate(Sample = "Sample 10a", CFA_CoughSneeze = NA),
sample10b = cbind(ega_scores_10b, data_addprefix(as.data.frame(predict(m7_all, newdata=data10b)), "CFA_")) |>
add_empirical(data10b) |>
mutate(Sample = "Sample 10b")
)
save(scores, file = "../data/scores.RData")